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Abstract 

£jj 1 We introduce a novel framework for developing statistical applications in health re- 

search, based on dynamic modeling of the investigated processes. We formulate the princi- 
ples of dynamic modeling in health research, which are coherent to those in other fields of 
research. Dynamic models explicitly describe causal relations which are to be adequately 
accounted in statistical methods, making them free of misuse of statistics and statistical 

■ fallacy. 
We propose the Dynamic Model of Population Health describing temporal changes in 

health indicators, having nature of state variables. The Dynamic Regression Method was 
developed as statistical method for the identification of the model. This method evaluates 
cohort trends for state variables at each age and calendar year. The method is illustrated 
by evaluating cohort trends for the Body Mass Index for men, using survey data collected 
in the years 1982, 1987, 1992, in North Karelia, Finland. 
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1 Introduction. 



Misuse of statistics and statistical fallacy are issues of concern in many fields, including medical 
research. The detailed classification of misuse and recommendations of how to avoid statistical 
fallacy could be found in books flJaffe and Spirer, 1987[ ), ( |Campbell, 19 74). One category of 



misuse, "lack of knowledge on subject matter" ( [Jaffe and Spirer, 1987[ ), could well be inter- 
preted as addressing causality among other things. 

Recently, it was acknowledged that a large proportion of published medical research contains 
statistical errors and shortcomings . "The problem is a serious one, as the inappropriate use 
of statistical analysis may lead to incorrect conclusions, artificial research results and a waste 
of valuable resources" ( Strasa k~et al., 2007] ) . Interestingly, the authors believe that one of the 
reason for misuse is lack of statistical literacy: "Medical researchers have to be encouraged to 
learn more about statistics, as various studies point to a lack of statistical knowledge among 
medical residents" and "statisticians should be involved early in study design" 
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In disagreement with this, we consider the case of statistical misuse, which occurs when 
a formally correct statistical method is used, however, causality is missing. These are the 
methods of evaluating secular trends in health indicators using data from a set of independent 
cross-sectional surveys (review on these methods is presented in section [5]). The linear secular 
trends are also used as a tool for intcrecensal estimates of population size (review of this 
could be found in ( Mo ltchanov et al., 1999[ ). We use term "secular trends" to refer to all these 
methods. 

Formally, methods evaluating secular trends look correct: data, linear model and assump- 
tions on data properties, in combination lead to evaluation of the model's parameters and 
testing hypotheses. In this schema, the wrong element is assumptions on data, usually sug- 
gesting smooth-line type of dependency of estimated means. 

Note that Hill's criteria for causation (see, for example, http:/ /en. wikipcdia.org/wiki/Bradford- 
HilLcriteria ), though sound reasonable, provide only circumstantial evidence for causality, and 
leave plenty of room for subjective judgments. 

To our view, the data on population size shown on Figure Q] provide a very strong evidence 
that secular trends do not exist in nature. Rather there are smooth evolution of population 
size along birth cohort. However, some of the statisticians and field researcher may disagree 
with this statement. 

To address properly and unambiguously causality in a real world process, first of all, knowl- 
edge and skill on dynamic modeling should be applied. Only at the next stage, statistical tools 
are to be considered to evaluate parameters of the dynamic model. 

Our practical task (target task) is to evaluate temporal changes in health on population level, 
using data from a set of independent cross-sectional population surveys. Traditionally, this task 
is approached by evaluating secular trends, which appeared to be a statistical fallacy. To build 
up an alternative approach, we develop Dynamic Model of Population Health (DMPH) and 
statistical method for its identification, the Dynamic Regression Method (DRM), producing 
time trends for health-related indicators within birth cohorts (Cohort trends, or C-trends for 
short). 

In turn, to build up a dynamic model, first we derive some principles, which we call the 
Principles of Dynamic Modeling in Health Research. These Principles are independent of the 
target task, so they could be applied to any other task, for example, to follow-up analysis with 
end-points. 

The aims of this paper are as follows: 

• to derive the Principles of Dynamic Modeling In Health Research 

• to develop the Dynamic Model of Population Health 

• to build up the Dynamic Regression Method and algorithm 

• to run the illustrative analysis 

Note that each of the four parts above is worth of more detailed, separate presentation. 
Therefore, the challenge was to provide concise and logically completed descriptions of all 
parts, clearly outlining logical interrelations between them. 

The earlier version of the Dynamic Regression Method was developed and presented by 



( Moltchanov and Mik'halskii, 2008 ), where C - trends were suggested as an alternative to 
circular trends, used so far. Historically, the method developed in MONICA for checking con- 
sistency of the reported demographic data ((Moltcha nov et al., 1999[ )) served as the prototype 



for the method developed by ( |Moltchanov an d Mik'halskii, 20081. Some general aspects, such 



as criteria for commonly used health indicators to serve as system State Variables, have been 
considered earlier by ([Moltchan ov, 1993 ). 
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Section [2] contains historical review of time-trend analysis. Section [3] presents history and 
key properties of Dynamic Modeling. Section [4] describes Dynamic Model paradigm being 
applied to health research on Individual and Population Level. In Section [5] we present the 
analytical model for the case of continuous, normally distributed one parameter. The example 
of application, employing the data of real study, is given in section [5J Section [7] contains 
conclusion and discussion. 



2 History of Time- Trend Analysis in Health Research 

Time-trend analysis of health indicators has been widely used so far in health research. One 
motivation for this comes from the fact that running a cross-sectional survey is the cost-effective 
way to collect the data for such an analysis (Mann, 2003). 

So far, this problem has been approached by assessing trends over time for means and other 
statistics (for example, percentiles) of the age-specific distributions of parameters of interest, 
such as traditional risk factor indicators (for example, systolic and diastolic blood pressure, 
cholesterol, body mass index), and their categorical derivatives (such as prevalence of high 
blood pressure, prevalence of high cholesterol, prevalence of obesity). 

Various terms are used in literature for such trends: "trends", "secular trends", "time 
trends". Here we will use the term "secular trends" for all of them. 

Among approaches used for trend analysis, the first one, "trends by linear regression", was 



the key element of the analysis in the WHO MONICA Project ( |Tunstall-Pedoe et al., 2003| ) 
Its steps include, first, calculating the age-group specific trends using linear regression, then 
aggregation using direct age standardization with fixed weights. This method has been applied 



to risk factors only (Dobson et al., 1998a| |Dobson et al., 1998b[ ), or to both, risk factors and 



rates ( |Kuulasmaa et al., 2000[ ). In the last case, the aggregated trends were subject to cor- 
relation analysis in order to test MONICA hypotheses. Some indication of problems in this 
approach come when the time plots of the age-group and survey-specific mean values of data 
items exhibit clear non-linearity of time plots and diversity of these plots over age groups (see, 
for example, ( |Tolonen et al., 2000| ), POL-TARa, BMI). 

A different modification uses the multiple logistic regression procedure applied to the whole 
set of data ( jGregg et al., 2005] ) . As a result, the marginal characteristics were obtained directly 
from the procedure. This is equivalent to direct standardization, with weights corresponding 
to the analyzed population. 

In examples above the method was applied to the samples having wide age range (40 years 
and more) while spans between consecutive surveys were 3-10 years. Some studies of adolescents 
deal with samples of age range 5-8 years, being sampled every year or every other year. In 
that case, the trends were first examined visually by age , since, as it was acknowledged, "they 
exhibit a wide diversity of age-specific patterns" ( |Kautiainen et al., 2002[ |Chen et al., 2003] ) . 

One commonly used approach to cope with such a diversity subdivides the overall time 
period into several segments and the overall age range into several age categories, for which 
the corresponding plots suggest linear trends. Alternatively, trends are evaluated for aggregated 
(age-standardized) parameters ( |Kautiainen et al., 2002| ). 

Summing up, we may conclude that methods used so far for the analysis of the changes 
of health-related indicators, though being the best available ones at that time, suffered from 
one principal drawback: lack of causality. In turn, this is due to the fact that comparison is 
made between different entities, or, in terms of Dynamic Model, between different objects (see 
section [3]) . 



3 



3 Dynamic Modeling: History and Inventory of General 
Principles 

In order to apply adequately the dynamic model paradigm to population health, we have to 
summarize the key notions and principles of this paradigm. We will start our inventory with 
the Newton's second law, which historically is the very first dynamic model. We will end up 
this section with formulation of the main principles of Dynamic Modeling in form of definitions, 
labeled PI - P4. 



3.1 Newton's Second Law of Motion 

The following stuff is quoting from the English translation by Motte ( |Newton and Motte., 1995] ) 
of the Isaac Newton's formulation in Latin in Principia, 1687. 

"Law II (in English): The alteration of motion is ever proportional to the motive force 
impress'd; and is made in the direction of the right line in which that force is impress 'd." 

We would like to use the structure of the above law's formulation as golden standard for 
general definition. In this view, the following comments are essential. 

• The above law was expressed verbally. The mathematical expression for it, for example 
in modern form 

is not fully equivalent to the original verbal form, since it does not specify what is cause 
and what is effect. 

• The law is not a kind of assumption to be used later in calculation, rather this is property 
of real world, derived by Newton from observation and logical inference. 

• The law postulates causality clearly: the cause is (motive) force, while the effect produced 
by this cause is the alteration of motion. 



3.2 Definition of Dynamic Model in General Case 

Mathematical forms of dynamic models, often being considered as a complex of several in- 
terlinked dynamic models, are subject to special mathematical discipline - theory of dynamic 
(dynamical)systems (see, for example, ( |Luenberger, 19 79)). It is important to stress, that our 
prime interest is in dynamic models, not in dynamic systems. 

Usually, there were no problems in application of dynamic models in engineering and econ- 
omy, as well as in interpretations of the results obtained. The need for re-inventory of the 
notion of dynamic model, has been encountered, however, when applying it in biology. 

The following definition was given by ( |Ellner and Guck cnhcim er, 2006[ ) : " Dynamic models 
are simplified representations of some real-world entity, in equations or computer code. They 
are intended to mimic some essential features of the study system while leaving out inessentials. 
The models are called dynamic because they describe how system properties change over time." 

We formulate the following definition: 

PI: Definition of Dynamic Models: Dynamic models arc simplified representations of pro- 
cess of change (over time) of some real-world entity, in verbal and mathematical form. 
They are intended to mimic some essential features of the study object in process of 
change, while leaving out inessentials. The models are called dynamic because they de- 
scribe changes of object properties over time, cased by driving forces. ( from Greek 
dynamikos " powerful" , dynamis " power" ) ; 



4 



Thus, dynamic model operates with object, characterized by set of properties, and Driving 
Force, acting at object and causing changes of the object's essential properties. 
Definition PI differs from one cited above in the following elements: 

• PI defines that dynamic model represents process of change of real- world entity, rather 
than real-world entity itself. 

• In PI "computer code" is excluded, since computer code may be an image of any process, 
possibly violating physical laws. While we are concerned with exploration of real world. 

• The models are "dynamic" not at all because of properties are changing. Rather because 
of the cause is postulated, generating these changes, - Driving Force. 

P2: Definition of Object: The real-world object is physical entity, traceable in time, which 
means that it could be observed over some time interval, possibly small enough, being 
virtually the same in common sense, and carrying the set of properties. 

P3: Definition of State Variables: An object is associated with a set of potentially mea- 
surable characteristics - State Variables. These characteristics are assumed to be essential 
properties of the object. For each such variable, the rate of change is proportional to the 
variable-specific net driving force, acting upon the object. Hence - State Variables are 
continuous and right-differentiable function of time. 

P4: Definition of Controls: The model ([IJ could be extended by adding position vector 
and written in "Dynamic Systems" format: 

dxi dx 2 1 , . , , 

-7T = x 2 , -T- = — u(t , 2) 

at at m 

where xi is position vector, x 2 - velocity vector, u(t) - represents driving force F , 
underlying the fact, that it could be non-continuous function of time. 

Note that in Driving Force for State Variable vector Xi is State Variable vector x 2 , 
thus being continuous, while Driving Force for State Variable vector x 2 is external force 
defined by (vector) function u(t). We will use term "Controls" to call such an external 
Driving Forces. 

For a body moving in gravitational field of a planet, State Variables position and velocity 
"predict" further behavior of this object, such as ability to "escape" the gravity of the planet 
without additional propulsion. In this respect State Variables differ from controls: the last ones 
modify State Variables rates of change, however they can not serve as predictors. 

Using dot notation for differentiation and combining vectors Xi and x 2 into common vector 
x, equations ([2]) could be rewritten as: 

x = ^(x,u(t)) (3) 

In verbal form, it postulates that rate of change of State Variables is a function of State 
Variables and controls. 

For any real- world object, a dynamic model is built-up to facilitate the following three main 
practical tasks pertaining to this object: 

Taskl: Analysis: Ultimately, design of function T is made by a researcher, dealing with a 
real- word object and its dynamics, depending on his/her intuition, experience and skill. 
The known prototypes play also important role. After design is made, the functional 
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form is known up to a set of still unknown parameters. The simplest case is linear form 
with coefficients to be set up. To evaluate these parameters, the measurements on State 
Variables and controls should be used, collected during some time period. We will use 
term "Analysis" for this task. 

Task2: Prognosis: If function T is identified and State Variables are known for a moment to 
and functions u(t) are defined on interval [t a , T] then State Variables could be calculated 
for this interval, thus making prognosis of future behavior of the object. 

Task3: Control: If criteria are set up, identifying the favorable future behavior of the object, 
and options are given for possible choice of controls as functions of time, then the task 
of control is to find out such a functions that optimize the given criteria. 

4 Dynamic Modeling in Health Research 

4.1 Health Research, Individual Level 

In the following example we consider a human-being whose weight change over time is subject 
of some study. Let x(t) be a weight of the subject measured in standard way at moment t. 
Due to measurements error it has a nature of random variable. In addition, the theoretical 
plot of this function over time will expose daily cycles, which are not of the study primary 
interest. Rather we would like to operate with some smoothed characteristics of weight. We 
assume that according to study protocol, weight is measured every day at the same time in 
the morning, after "emptying body tanks" before eating. So, we may think of sequence of 
measurement time moments f,-, where i is sequence number of day since the beginning of the 
study. 

We introduce function v(t) defined as follows: 

v(U) = E(x(U)), t = U, 

v(t) = v(ti) + (v(t i+1 ) - v(U)) ■ ~ 4 U<t<t l+1 (4) 

H+l — ti 

This function satisfies all the requirements for the State Variable: it is continuous and right- 
differentiable function of time. The current knowledge on weight changes in adults could be 
summarized as follows. Weight change in adult is, in fact, change in amount of body fat, which 
is determined by balance of calories taken with meal and burned throughout the body activity 
in over a certain time period. Thus, we can introduce function u(t) representing daily balance of 
calories, expressed in weight units (see, for example, http: / /www. weightlossforall.com/calories-per-pound. htm 
"One pound of body fat equals roughly 3,500 calories."). This function plays role of control 
for v(t) and we may postulate simple model for weight change: 

v = u{t) (5) 

We consider a hypothetical study testing some technique for weight reduction (it may include 
education, dietary recommendation, advice on physical activity etc.). Assume that measure- 
ments of weight are available before and after the beginning of intervention, moment to- 

To highlight principal conceptual aspects, we make the following additional simplifying 
assumption: u(t) — u±, if t < to, u(t) — U2, if t > to, where u\, ui are scalars. In practice, 
estimates for U\ and u% could be obtained as slopes in linear regression models applied to 
measurements x(t) for t < to and t > to correspondingly. 
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Condition u-i > u\ indicates that tested technique is better than previous one (in practice, 
value a could be added having sense of " practical significance" , so that condition will look like 
U2 > tii + ot). Note, that "success" is derived from comparison of time trends for weight, not 
from the fact of decreasing weight for t > to- Theoretically, positive u-i may indicate success, 
if the weight growth has diminished, and negative ui is not a success, if u\ also was negative 
and approximately the same in value. 

It is convenient to call all the data items, available in the study database and pertaining to 
a subject at certain moment of time, measurements. 

In dynamic model view, most of the measurements fall into three categories: 

1: State Variables: for example, age, weight, height, schooling years. Recall, that mea- 
surements for State Variables are not State Variables. They refer to each other as x(t) 
and v(t) in the described above example. 

2: Modifiers: for example, smoking status (smoking now, ex-smoker, or never-smoker) , cur- 
rent physical activity, current dietary habits ( including 24 hours food consumption 
recall) 

3: Class indicators: for example, sex, race, community, other characteristics, which are 
categorical and believed to be constant during the study time span. 

Outside of the above categories are multi-item outcomes of different questionnaires and 
tests. Some of them could result in one summary item ( for example, current physical activity 
level, which then could be classified as control). Questionnaire on smoking history may result 
in total amount of tobacco smoked so far. This indicator has a nature of State Variable. We 
leave further consideration of this issue for future publications. 

Similar to models in mechanics, Modifiers in health research modify status of body in terms 
of State Variables, however they can not serve as predictors, for example, as predictor of 
instantaneous failure. Thus, current Hazard can not depend, for example, on smoking or phys- 
ical activity. This simple rule of dynamic model philosophy is widely violated in practice of 
methods used so far in health research. 

Observe that continuous function of any number of State Variables is itself a State Variable, 
and the same function applied to the measurements of corresponding variables serve as mea- 
surement for resulting variable. The expression ([S]) remains valid for this variable after control 
u{t) is properly scaled. 

In our future example we will deal with such a variable, The Body Mass Index, defined as 



BMI= ^™ (6) 



weight(kg) 
height(m) 2 

Here we pay tribute to tradition, using term weight instead of scholastically more correct term 



4.2 Dynamic Model of Population: Heuristic Approach 

We may think of population as a collection of subjects identified for each calendar date by 
a certain rule. For example, for population of an urban district such a rule may identify all 
permanent citizens having home address within unambiguously defined administrative bound- 
aries. The rule must be the same throughout the calendar period for which the population is 
supposed to be analyzed. 

Health-related and other population characteristics, if available, has a form of age-distributed 
profiles, specific for calendar years, rather than individual-specific measurements for all current 
subjects of population. 
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Measurements on population level are performed for random samples (stratified or not), 
taken, for example, every 5 year. 

Thus, the challenge is how to adopt for population level the dynamic model paradigm 
described so far for individual level. 

To describe population history, it is convenient to use plane (y,a), where y is real- valued 
calendar time in years, vertical axis, a is real-valued age in years, horizontal axis. Such a set 
up for axes anticipates further use of matrices with indexes y, a, when the first index is row 
number (vertical coordinate). For consistent setup, we have to specify an observational frame 
in terms of ranges \y m in,Vmax] for V and [a min ,a max ] for a. 

Each subject may enter this population due to birth (if a m i n = 0), or crossing left-low 
boundaries, or migration in. Each subject may leave this population due to death, migration 
out or crossing the right-upper boundaries. If a subject with coordinates (yo, ao) is within the 
population during time t , at that time it has coordinates (yo + t, ao + t). Thus we may say, it 
is moving along cohort line. 

Consider all subjects having coordinates on half-open interval ({yo, ao — Aa), (yo, ao)] at time 
t = 0. At time At all those left in population will arrive at ((yo + At, ao — Aa + At), (yo + 
At — Aa, ao + At]. In other words, the birth cohort of width Aa moves from (yo,ao) to 
(yo + At, ao + At). We may think of such a cohort as of a container moving on plane (y,a). 
The contents of each container in process of movement is changed due to migration and death. 
If the rate of contents update is negligible (say, less than 1% per year), we may ignore it in 
our analysis. If not, the analysis has to take this into account. 

Each container fits the definition of the dynamic model object, if we regard the correspond- 
ing State Variable as mean of State Variables for currently available subjects. The dynamic 
equation then could be obtained from ones for each subject, having form |5]>. by taking means 
of both sides: 

ij=u(t) (7) 

Since the whole selected observational frame could be covered by collection of non-overlapping 
cohorts of selected width, we may conclude that, in case of population, the overall dynamic 
model is a collection of dynamic models specific for each cohort. 

4.3 Dynamic Model of Population: Axiomatic Setup 

The theoretical abstraction for birth cohort is one of infinitesimal age range, characterized by 
multidimensional distribution of the parameters of interest, not by physical subjects. 
Let C be a 2-dimensional real compact: 

C {(y^a) . y \ymin,ymax\i ® ^ [^mm^maz]}) 

where y is calendar time in years and a is age in years. 

Consider a population defined on this compact, which suggests that there potentially exists 
a set of random variables (r.v.) X^, i = l...k representing the corresponding set of measurable 
indicators of interest (State Variables) defined at each point (y, a) of compact C. In this paper 
we restrict ourselves to the case of one indicator, so that subscript of X will be omitted. To 
make the following description more illustrative let us keep in mind the Body Mass Index 
(BMI) as an example of the indicator in question. 

We introduce the following notation 

v(y,a)=E(X(y,a)). 
For the sake of simplicity while describing the core dynamic model, we assume, 
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X(y, a) = v(y, a) + e, where E(e) = 0, D(e) = a 2 , V(y, a) : (y, a) e C (8) 

The dynamic equations describe changes of the distribution of r.v. X for a birth cohort 
taken at point (y, a) over time interval dt: 

v(y + dt, a + dt) = v(y, a) + u(y,a)dt + o(dt), where ^— ^ -)• 0,asdt 0. (9) 

dt 

On one hand, function u(y,a) represents the rate and direction of change of the State 
Variable due to the driving force generated by the environment. On the other hand, it is the 
driving force (control) itself, properly scaled. 

The driving force at {y,a) does not depend on the properties of the cohort passing at the 
time y the age a. Moreover, theoretically, the very fact of its existence doesn't depend on 
whether or not there is a non-empty cohort passing at the time y the age a. 

For the sake of convenience we will use terms "Mean levels" or "levels" for the values of 
function v(y,a), and "cohort trends" or "C-trends" for the values of function u(y,a). 

In the advanced model the function u(y, a) represents sum of the environmental force and 
the force due to current state of the cohort. This will lead to replacement of u(y, a) in © by 
u(y, a) + bv(y, a), where b is a model parameter. 

Let vo(y, a) be the value of v(y, a) at low-left boundary of the compact C for a (birth) cohort 
crossing the point (y, a): 

v (y,a) = v(y - S, a - 5), where S = min(y - y min , a - a min ). (10) 

Then v(y, a) can be expressed as 

v(y,a) = v (y,a) + u(y-t,a-t)dt 
Jo 

Thus, if the values of vo(y,a) at low-left boundary and u(y,a) on C are known, then the 
function v(y,a) could be evaluated for each point on C. 

The generalization of the model (JSJ), © for the case of multidimensional distribution and 
state-dependent dynamics is straightforward, by treating functions v(y, a) and u(y, a) as vector 
functions, by replacing D(e) — a 2 in (JSJ) by Cov(e) — £ and by replacement of u(y,a) in © 
by u(y, a) + bv(y, a), treating b as a matrix. 



5 Dynamic Model of Population: Analytical Form 

5.1 General Formulation of the Task 

Suppose that a set of measurements is available (xk, yt, Qfe), k = 1, . . . , K , for subjects selected 
in a set of the independent cross-sectional surveys. We assume that for each survey the stratified 
by gender and age group random sample scheme was used. The age group stratification could 
be different in different surveys, however, for standard case, we assume that overall age range 
is the same in all surveys. 

The general formulation of the task is to estimate the functions wo(y,a) and u(y,a) on C, 
using the available measurements {xk,yk, o-k), k = 1, . . . , K . 

To solve this problem one option would be to formulate the optimization problem in func- 
tional space: to minimize the functional /: 
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I(u,v )= (^2(x k -v(y k ,a k j^J 



(11) 



applying some additional requirements on functions u(.,.) and i>o(-,.), such as continuity 
(piece- wise continuity), and /or restricted variation. 

However, it seems more convenient to transform the above problem into the discrete - scale 
analogue and to take the advantage of the simplicity of the analysis and adaptation of the 
numerical methods available in the standard statistical packages. 

5.2 Discrete-Scale Model 

Let i and j be an integer value of time in years and an integer value of age in years corre- 
spondingly. Our intention is to build up the integer- values proxies of the equations ([8] - ITT|) . 
Let P(i,j) be a parallelogram-shaped element (convex hull) defined by its angle points: 

{(ij-i), (ij), (i+u+i), (i+i,m 

excluding its left and upper boundaries, which could be written as 

P(i,j) = {(a,y):ye[i,i + 1), a e ((j - 1) + (y - i), j + (y-i)]} (12) 

We impose for function it(.,.) the conditions of being constant on each P(i,j) and for 
functions «(., .) being constant on a and linear on y with constant slope u(i,j). 
Formally this could be expressed as follows: 

u(y,a) =u(i,j), Vi,j,y,a: (y,a) e P(i,j) (13) 

v(y,a) = u(i,j)-(y-i)+v(i,j), Vi,j,y,a: (y,a) € P(i,j) (14) 

We derive minimal and maximal values for i and j from the correspondent values for y and 
a using definition (IT^I) : 

(^mmjimrn) • (?/mm?^mm) £ P{^min, jmin^) , (1^) 
(imax , jmax*) • {]Jmax-} Q>max) £ P{^raax, jmax) 

For convenience, from now on we will use relative scale for age and time, defined by trans- 
formation 

^ ^min ^ ^) J Jmin ^ J 

Consider functions u(i,j) and v(i,j) defined on integer- valued two dimensional domains 
U = [0,1], j€[0,J]}, 

V = [0,1+1], je[0,J + l]}, (16) 

correspondingly, where 

I imax irnin, J jmax jmin 

Now the main dynamic equation ((9] ) could be rewritten as 

v(i + l,j + l) = v(i,j)+u(i,j), V(i,j)€U (17) 
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Let vo(i,j) be the value of v(., .) at low-left boundary of the domain V corresponding to a 
(birth) cohort crossing the point 

v o{hj) — v(i — J — S), where S — min(ij). (18) 

Combining (|17p and (fT8|) . we rewrite equation (jlOl) as: 

8 

v(i,j) = v (i,j) +^2u(i- m,j - m) (19) 

m— 1 

From (fTO)) it follows that if v(i,j) is set up on the low-left boundary of V and u(i,j) is set 
up on the whole W then v(i,j) could be calculated for the whole V. 

Finally, assembling (|5j). IT9"| and (fLT| for each available observation (xk,Uk, CLk), k = 1, . . . , K, 
we obtain: 

<5 

Xk = v (i,j) +^2u(i- m,j - m) + (y k - i) ■ u(i,j) + e k , 

m—1 

where Var(e k ) = a 2 , Cov(e k , e{) = 0, if k / 1 (20) 
Let z be a vector with components vo(i,j) and u(i,j) ordered in the following way: 

v = (v(I + l,0),...,v(0,0),...,v(0,J + l)f 
u = (ii(0,0),...,u(0, J),...,u(I,Q),...,u(I, J)f 
z = (vj I " T ) T (21) 
Using vector z and introducing vector of coefficients b^, we can rewrite (|20l) in the form 

x k = (b fc , z) + e k , where Var(e k ) = cr 2 , Cov(e k , ei) = 0, if k ^ 1 (22) 

This form represents a particular case of Gauss-Markov Setup for the Least Squares Linear 
Estimation problem ( |Rao, 1973] ) . 

Let Bo be a matrix composed of row vectors in (|22[) . z and Xo stand for column vectors 
of the parameters Zj and the variables Xk correspondingly, and So be a scalar function defined 
as 

Sb(z) = (B z-xo) T (B z-xo) 

Note that if ranfc(Bo) = rfim(z), then estimates obtained by unconditional minimizing of 
function <So(z) are unique ones. Such a case takes place only if the observations cover all the 
elements P(i,j) when surveys cover the whole analysis period without gaps. 

In practical cases, minimizing of So results in singular or ill-posed Inverse Problem, and 
so-called regularization techniques are needed to obtain meaningful solution estimates. Most 
of these techniques employ the idea of smoothing of some function having clear physical inter- 
pretation ( |Neumaier, 1 999). 

Here we suggest one such technique for smoothing. 
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5.3 Smoothing 

We define the following indicator of smoothness of function v{., .) 



Si(z) = El=o E/=i («(i,J - 1) - 2«(i, i) + i + 1) 
E/=o Eli Ui - - 2v(i,j) + v(i + l,j) 



(23) 



Each term in this sum represents the square for a proxy of the second derivative of function 
v(., .) with respect to age or with respect to calendar time at point 

Replacing v(.,.) by vq(.,.) and u(.,.) using fT^]) , and the last ones by vector z, we will 
transform the previous expression to the following form: 



Si (a) = (B lZ -0) T (B lZ -0) 
Similarly, we define indicator of smoothness of function u(., .) 

S 2 (z) = Ej=oE/=i - 1) - MU) + u(i,j + I))' + 



(24) 



allowing form 



E/=o E-=i («(* - U) - 2«(i,i) + «(< + 1, j) 



S 2 (z) = (B 2 z-0) T (B 2 z-0) 



(25) 



Now we can add one or both constraints Sfc(z) < Qfc with some selected a& > 0, k = 1, 2, to 
the model (|22p. Observe that indicators So, Si, S2 are quadratic functions in finite vector space 
E„ with elements (vectors) z and n = dim(z). The optimization problem for point estimation 
for our case, could be formulated as 



min S*o(x), subject to Sfc(x) < au, with given Oik > 0, k = 1, 2. 

x6E„ 



(26) 



Let no , ni and ri2 be numbers of rows in matrices Bo , Bi and B2 correspondingly. Let Ai, 
A2 be some non-negative scalars. Introducing matrices and vectors 



B 



B 



B, 

B 2 / 





V J 



w 



Io 











A1I1 











A2I2 ( 



(27) 



where Io, Ii and I2 are identity matrices of rank uq , n\ and ni correspondingly, we can 
formulate the problem of least squares estimation in the following form (a modification of 
Gauss-Markov setup which fits form of Aitken setup ( |Rao, 1973] ) 

x = Bz + e, E(e)=0, D(e) = ^W" 1 (28) 

for which the point estimation problem is 

min S(z), whereS(z) = (Bz - x) T W(Bz - x) = S (z) + AiSi(z) + A 2 S , 2 (z) (29) 

(Moltchanov and Mik'halskii, 2008) have shown that problems (|2"o| and (12^1) are equivalent: 
problem (|2l)]) with given a.\, possesses the same solution as problem (|2T)]) with some Ai, A2, 
and vice versa, or both don't possess any solution. 
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Since part of its components are set to zero, the data vector x in (071) could not be treated 
as a "true" data vector if the problem is considered from the classical frequentist prospective. 
The last one is based on retrospective evaluation of the procedure used to estimate parameters 
over the distribution of possible data values conditional on the true unknown values of param- 
eters ( |Gehnan et al., 19 95), p. 7. The logically consistent treatment of the problem is based on 
Bayesian paradigm, where statistical conclusions about unknown parameters are made in terms 
of probability statements, conditional on observed data. As noted in ( Gelman et al., 19 95), p. 7, 
in despite this difference , it will be seen that in many simple analyses, superficially similar 
conclusions result from the two approaches to statistical inference. In particular, this con- 
cerns Bayesian analysis of the classical regression model: under a standard noninformative 
prior distribution, the Bayesian estimates and standard errors coincide with the classical re- 
sults (Gelm an et al., 19 95), p. 235. The last statement justifies use of classical formulas and 
numerical procedures for our "non-classical" case. 

The question of primary practical importance is the existence of a unique solution for the 
problem ([29)), 

The following statement is proofed in ( |Moltchanov and Mik'halskii, 2008] ) : 



Corollary 1 For existence of a unique solution to problem (l29l) it is sufficient to have 4 data 
points such that the corresponding points (y, a) on plane y, a satisfy condition: no any 3 of 
them are located on a common straight line. 



5.4 Outlines of the Algorithms. Setting up the Regular izat ion Parameters 

As soon as parameters Ai, A2 are given in setup (f2"Tl |2"5|) . the following could be obtained 
routinely: z - point estimate of vector z, covariance matrix of this estimate Ccw(z), and a 2 - 
estimate of a 2 . 

Using functions u(i,j) and v(i,j) defined in l|16p . we can create matrices 

V:v itj =v(i + l,j + l), 

U : Uij = u(i + 1, j + 1), 

and vectors 

v = (Shape{V, 1)) T , 

u = (Shape(U, 1)) T , 

where Shape is matrix function reshaping the original matrix into resulting one with different 
number of rows and columns (available, for example, in SAS/IML ( |SAS Institute Inc., 2004b| ). 
In our case, results are vectors with consequently concatenated rows of the original matrices. 

Each element Vij of matrix V corresponds to element Vk of vector v with 

k = (i - 1) ■ ncoZ(V) + j, (30) 

where ncol(.) is matrix function returning number of columns. Similar rule could be applied 
for linking U and u. 

Definition of vector u in (f2~T|) and expression for functions u(i,j) in (TT9"]) allow to construct 
matrices 

A z2v : v = A z2v z, and 
A z2u : u = A z2u z. 

Hence, the covariance matrices could be derived as 
Cov(u) = A z2u Coi;(z)A z2u T , 
Ccw(v) = A z2v C'cw(z)A Z 2v T , 

from which the corresponding matrices of Pearson's correlation coefficients R v and R u could 
be routinely produced. 
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Consider two consecutive level estimates v(i,j), v(i,j + l) allocated along age axe (similar 
consideration could be applied to allocation along calendar years. 

The coefficient of correlation for these estimates could be derived from R v applying rule (|50)l . 
Let denote it r a ,i,j- Similarly, coefficients of correlation r Vt ij could be derived for estimates 
v(i,j), v(i + 1, j) allocated along years axe. 

Consider task of predicting estimate v(i,j + 1) using linear predictor based on v(i,j). The 
expressio n 1 — i is proportion of "unexplained" part of variance of v(i,j + 1), (see, for 
example, ( |Rao, 1973] ) p. 266). This part could be interpreted as "new information", or "signal", 
while r 2 i ■ could be regarded as proportion of " Noise" . The better smoothness is associated 
with lower signal-to-noise ratio. We combine all local indicators of smoothness into one common 
vector 

f v = Shape(V smSL , l)\\Shape(V smy , 1), where v sma ,i,j = 1 - j, v smytitj = 1 - r* tiij 

(31) 

Vector f u could be defined in similar way. 

For practical use we have to select function, producing sample statistics for a vector- 
argument, such as mean, median, minimum or a value of one predefined component and a 
target value for this statistics, f sm - Let fstat be generic name for such a function. Then 
iterations are run by selecting Ai and A2 until the following condition is satisfied 

max(abs(log(fstat v (f v )) - log(f smv )),abs(log(fstat u (f u )) - log(f smu )) < 5, (32) 

where 5 is predefined accuracy. 

With increasing values of f sm v, fsmu the corresponding lines and surfaces visually become 
smoother. For level estimates, for example, if f sm v ^ 0, then Ai — > 00, and solution converges 
to 4-parametric surface (( |Moltchanov and M ik'halskii, 2008] )). 

To measure difference in C-trends over age and calendar year, the pairwise comparison tests 
are performed for mean values of C-trends, evaluated for a set of age-year clusters, defined by 
cluster sizes, A Q and A y . 

Let U c be matrix of such mean values, u c —Shape(XJ c , 1) T and matrix A U 2 U c such that 
u c = A u2 ucU. As soon as, matrix A u2u c is created for given A a , A y , U c could be calculated, 
as well as variance/covariance values for its elements in a format of 

C = Cov(u c ) = A U 2ucCcw(u)A U 2uc T - 

For each cluster, statistics and corresponding probabilities are computed for pairwise com- 
parison of mean C-trends for current cluster and for adjacent one for older age group, and 
for current one and for adjacent one for the next calendar years period ( if the corresponding 
clusters exist). Using classical paradigm, this is done by testing linear hypotheses in form 

H : u ci - u cj = 0. 

General expression for F- value (see for, example, SAS/Stat manual, QSAS Institute Inc., 2004c| )) 
in this case takes a simple form 

F _ {Uci-U CJ ) 2 



Corresponding probability is computed using SAS function probF (see QSAS Institute Inc., 20lT] )) 

as 

Pr = 1 — probF(F, 1, n — r) 

Note, that in Bayesian view, these probabilities should be referred to as tail-area probabilities 
for posterior predictive distributions (( |Gelman et al., 19 95) , p. 169). 

The algorithm, implementing the above outlines, is written in SAS code using SAS products 
( dSAS Institute Inc., 2011] ), citeSAS9.2PROC, ( |SAS Institute Inc., 2004b] , ( |SAS Institute Inc., 2004c| ) 
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QSAS Institute Inc., 2004a[ )). Results of pairwise tests are presented graphically in figure, pro- 
duced by PROC GCONTOUR, properly annotated ( see Figure [5] in example of application). 

For reference, we will call this algorithm DRM2(R), with prefix DRM2 to differentiate it from 
those developed in ( Moltchanov and Mik'halskii, 2008 ). We have built up also the modification 
of this algorithm, processing aggregated data, DRM2(A), thus DRM2(R) for "oRiginal" data, 
and DRM2(A) - for "Aggregated" data. 

Original individual data may be of quite big size, which reflects row number no of matrix 
Bo in (1271) . and hence, the required memory and time for calculation. 

Aggregation is applied to original measurements (xk,yk, Ofc), k = 1, . . . , K , producing sum- 
mary statistics for {age ■ year) cells with size 1. As a result, arithmetic means are produced 
(x c ,y c ), number of original measurements in each cell (n c ), and Scss c -Corrected Sum of 
Squares, where c — 1, . . . , C - collection of non-empty cells. Matrix Bo and vector Xo in (|27l ) 
should be replaced by Bo and Xo, with cell-specific rows. 

Let no be a frequency vector with components (n c ). The following expressions are essential 
elements of the DRM2(R) algorithm. 

Contribution to cross-products B ■ x and B ■ B: 

EXPR1: (B T • Diag(n ) ■ x ) 

EXPR2: (B T • Diag(n ) ■ x ) 

Contribution to sum of squares of error terms, 

EXPR3 : (B z - x ) T • Diag(n ) • (B z - x ) + £f=i S C ss c 

Note, that DRM2(A) and DRM2(R) will produce identical outputs if all within cells are 
equal. 



6 Example of Application 

6.1 Data 

To illustrate the method and to demonstrate its performance, the data from three cross- 
sectional surveys, conducted in North Karelia, Finland, during the period 1982 -1992, will be 
used. Formally this set of data can be characterized as follows: 

• Study population: North Karelia, Finland. 

• Study period: 1982-1992; 

• Source of data: cross sectional independent surveys conducted in years 1982, 1987, 1992 

• Sampling frame for each survey: the stratified by 10- year age groups (25-34, 35-44, 45-54 
and 55-64) and gender random sample. 

The following specifications defines sub-sample of records and items selected for analysis. 

Only data for men will be used; the number of examined men in years 1982, 1987, 1992 is 
equal to 1537, 1481 and 673, correspondingly. 

Original measurements of interest are: gender, date of birth, date of examination, weight 
and height. 

The analysis variables included in the model: 

BMI - the Body Mass Index, defined as weight(kg)/height(m) 2 . 

AGE - age in full years, defined as year of examination minus year of birth. 

YEAR - date of examination measured in years. 

All surveys have started at the beginning of the year, surveys 1982 and 1987 have been 
completed in 4 months, survey 1992 - in three months period. 
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6.2 Analysis Setup 

The algorithm modification DRM2(A) have been used, preprocessing original data into aggre- 
gated format. There were 120 aggregated observations, 40 for each survey year. The analysis 
was set up for the age range 25-64 and for the calendar year period 1982-1992. In rule, con- 
trolling iterations, (1321) . smoothing factors f sm v, fsmu were set to 0.2, accuracy level, S was 
set to 0.05; cluster sizes, A a , A y , for producing comparison tests were set to 5. 

We have found that for practical purposes it is enough to use one selected point, (1,1) for 
fstat v , and one selected point, (int(ncol(\J)/2),int{nrow(\J)/2)), for fstat u . 

6.3 Results 

The results of the analysis are visualized by the set of 3-dimensional figures. 

Figure displays the values representing means of BMI calculated for each age and year, 
for which the survey data are available (number of cases in each cell exceeds 9). To visualize 
the along-cohort changes, the columns corresponding to the same birth cohorts in different 
surveys have similar shades of grey. 

Figure [3] displays estimates for the mean levels of BMI for the whole domain, with study 
age range plus one year, and study period plus one year. 

Figure [^displays C-trends with 95% confidence intervals, shown at left and front boundaries 
only. 

Figure [5] displays mean levels of C-trends for specified age- year clusters, with P-values for 
differences between clusters. 

These figures illustrate the principle "one figure is better than one hundred tables" , though 
all the underlying data are available and could be presented in a set of tables. 

Figure [3] shows that mean BMI levels increase along cohort lines throughout the study 
period, although they are different for different birth cohort. Specific peaks and troughs follow 
cohort lines. 

Recall that C-trends represent the net external Driving Force (Modifier) causing changes 
over time in cohorts. Therefore, changes in C-trends pattern over calendar years may indicate 
effect of preventive activities, while difference across age range may indicate both, age-specific 
uncontrolled changes and/or different susceptibility to prevention. 

In our case, Figure [5] shows clear decrease of C-trends in the period 1987-1992 compared 
with the period 1982-1986 in the age range 35-40 (p < 0.05); No other significant differences 
between adjacent clusters were detected. 

The further detailed analysis and final interpretation of the results may require a log of the 
events affecting the socio-economic and health care profiles of the study population during 
the study period. For example, a feasible explanation of the observed effect in C-trends in age 
range 35-40 could be associated with creating new working places in years 1987-1992, which 
have decreased population flow out of the area, taking place in years 1982-1986 in this age 
range and leading to negative health selection (subjects with low BMI were leaving the area 
in searching for job places). 

Summing up, we can conclude that, in general, clustering of C-trends looks reasonable, so 
we can use the results of comparing C-trends levels in adjacent clusters for our analysis. 

7 Conclusion and Discussion 

In this paper we have presented a novel formulation of the key principles of dynamic modeling in 
general, and in application to health research, which justify the structure and interpretation of 
the core models dealing with C-trends. In particular, according to these principles, traditional 
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risk factors' indicators fall into two categories, State Variables and Modifiers (see section [3]), 
having different dynamical nature and, hence, playing different roles in the model and analysis. 

As corollary of this, circular trends for State Variables have no sense at all. At the same 
time, only State Variables may determine instantaneous hazard rate of failure. In dynamic 
models, causality is postulated: changes are due to Driving Forces (Modifiers), existing in the 
real world. In case of consecutive survey data, C-trends are believed to be proxies for Driving 
Forces, providing the tool for three main practical tasks: analysis, prediction and control of 
health on population level ( see section [3]) 

We have used these principles as a framework for developing the dynamic model of simulating 
the temporal changes in characteristics of a real-world object - population. In the course of 
this process, first, we have identified two interacting objects, population and its environment, 
on the top aggregation level. Further system analysis has led us to breaking down the study 
population into a set of potentially infinitesimally narrow birth cohorts, carrying over time 
health state profiles expressed in terms of health related indicators (State Variables). 

The model employs the health field concept, suggesting existence of an influencing factors 
(Modifiers), generated by environment and acting on the population, specific for each calendar 
year and age, and causing within-cohort changes of the health indicator with rate of change 
corresponding to the strength of this factors. 

For illustrative purposes we have selected one-parameter case with continuous, normally 
distributed parameter and with strength numerically equal to rate of change. While keeping 
model reasonably realistic, these simplifications help to highlight the key properties of the 
dynamic model of population health and method of its identification - the Dynamic Regression 
Method. 

In the illustrative example, we have shown that the Dynamic Regression Method provides 
a sensible view on the BMI dynamics. It reveals clear difference between the levels of the 
parameter and its C-trends. From practical prospectives, it is C-trends, not levels, which 
primarily seem to be modifiable by preventive activities or involuntary changes affecting the 
population. It is worth noting that outcomes from the DRM analysis serve as data for the 
next-level analysis, involving other information and aiming at finding reasonable explanation 
of the observed dynamics (diagnostic property of DRM) . One of the important complementary 
component for such an analysis is dynamics of the population size ( we have developed a 
modification of the DRM for that type of data, this is a subject for one of the next publication). 
If there is significant migration "in" or "out" of the study population, the observed effects could 
be entirely or partially due to the population instability (health selective effect). The outcomes 
from the DRM analysis could be used straightforwardly for prediction of the age-specific profile 
of the State Variable, say, for 5 year period, by applying the C-trends at the last year of the 
study period to the estimates of the parameter's levels at that year. Such a projection will not 
cover the cohorts, not included in the study age range at the last study year. 

Recall that this method has been developed as an alternative to the secular trends used so 
far. In this respect, it is worth noting that the model presented here is characterized by local 
cohort trends (C-trends), which have clear interpretation: changes in the State Variable of the 
same physical entity per time unit. If we will formally calculate a characteristics resembling 
age-specific secular trend, we will obtain a difference between two different physical entities 
(birth cohorts), caught occasionally at the moments of measurement. Hence, it may behave 
quite arbitrarily. In other words, in the view of the dynamic modeling approach, secular trends 
do not exist in nature. In one special case only, when all the age profiles of a State Variable 
are the same over calendar years (stationary case), formally calculated secular trends will be 
equal to zero at each age within the study age range. Only in that trivial case, secular trends 
possess both, predictive and diagnostic power. However, even in this case, secular trends are 
kind of statistical fallacy, since missing causality. 
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There are certain restrictions in using the current version of DRM methods, imposed by the 
size of the problem, due to using matrix operations. Transfer to the Bayesian Data Analysis 
and using Markov chain Monte Carlo simulation methods (Gelman et al., 1995) seems to be 
a solution for these problems. 

The simplified dynamic equation used in the current model could be modified, accounting 
for the fact that rate of change may depend also on the current level of the State Variable. 

Finally, the most comprehensive model needs to be developed, comprising multiple State 
Variables, and corresponding C-trends as a linear functions of current State Variables. Such a 
model could be a powerful practical tool for prediction of population health for about 5 year 
span. 
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Figure 1: Population size, Men, Original data. Count by year and age. 
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Figure 2: BMI, Men, Survey data. Means by year and age. 
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Figure 3: BMI, Men, Estimates of means by year and age. 
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Figure 4: BMI, Men, Estimates of C-trends by year and age. 
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Figure 5: BMI, Comparison of C-trends by clusters of age and calendar years. 



